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Solid  oxide  fuel  cells  (SOFCs)  are  considered  to  be  among  the  most  important  fuel  cells.  However,  SOFCs 
present  a  challenging  control  problem  owing  to  their  slow  dynamics,  nonlinearity,  and  tight  operating 
constraints.  In  this  paper,  we  propose  a  model  predictive  control  (MPC)  strategy  based  on  genetic  opti¬ 
mization  to  solve  the  SOFC  control  problem.  First,  a  support  vector  machine  (SVM)  model  is  identified 
to  approximate  the  behavior  of  the  SOFC  system,  then  a  specially  designed  genetic  algorithm  (GA)  is 
employed  to  solve  the  resulting  constrained  nonlinear  predictive  control  problem.  A  terminal  cost  is 
incorporated  into  the  standard  performance  index  to  further  enhance  the  control  performance.  More¬ 
over,  the  GA  is  accelerated  by  improving  the  initial  population  based  on  the  optimal  control  sequence 
obtained  for  the  previous  sampling  period  and  a  local  controller.  In  addition,  a  dynamic  constraint  is 
also  adopted  in  order  to  meet  the  requirements  for  the  desired  fuel  utilization  and  control  constraints. 
The  measures  to  achieve  offset-free  properties  are  also  discussed.  Simulation  results  on  an  SOFC  system 
illustrate  that  the  proposed  method  can  successfully  deal  with  the  control  and  control  move  constraints, 
and  that  a  satisfactory  closed-loop  performance  can  be  achieved. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Solid  oxide  fuel  cells  (SOFCs)  are  considered  to  be  among  the 
most  important  fuel  cells.  Among  the  various  types  of  fuel  cell,  solid 
oxide  fuel  cells  (SOFCs)  have  attracted  considerable  interest  as  they 
offer  a  wide  range  of  applications,  flexibility  in  the  choice  of  fuel, 
high  system  efficiency,  etc.  [1  ].  However,  SOFCs  present  a  challeng¬ 
ing  control  problem  owing  to  their  slow  dynamics,  nonlinearity, 
and  tight  operating  constraints  [2]. 

Model  predictive  control  ( MPC)  appears  to  be  the  most  appropri¬ 
ate  control  strategy  for  SOFCs.  A  comparable  H ^  control  strategy  for 
an  SOFC  model  has  been  shown  to  be  unsatisfactory  [2].  Recently,  a 
data-driven  linear  MPC  strategy  using  subspace  system  identifica¬ 
tion  was  proposed  in  Ref.  [3]  to  control  an  SOFC  system.  However, 
more  researchers  have  resorted  to  nonlinear  MPC  strategies  for 
better  control  performance  [4-9].  These  studies  have  typically 
employed  nonlinear  models,  such  as  the  Hammerstein  model  or  the 
radial  basis  function  (RBF)  neural  network  model,  as  predictors,  and 
then  a  non-deterministic  polynomial-time  hard  (NP-hard)  nonlin¬ 
ear  optimization  problem  was  solved  on-line  using  conventional 
nonlinear  optimization  techniques,  e.g.  sequential  quadratic  pro¬ 
gramming  (SQP).  In  contrast  to  these  studies,  genetic  algorithm 
(GA)-based  MPC  has  been  widely  studied  in  recent  years  because 
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GAs  have  a  number  of  advantages  over  conventional  nonlinear 
optimization  techniques  in  solving  the  constrained  nonlinear  opti¬ 
mization  problem  [10-15]. 

In  this  paper,  we  propose  a  constrained  MPC  strategy  to  solve 
the  SOFC  control  problem  based  on  a  support  vector  machine  (SVM) 
and  genetic  optimization.  A  terminal  cost  is  incorporated  into  the 
standard  performance  index  to  make  it  possible  to  adopt  short 
horizons  while  maintaining  a  satisfactory  performance.  Moreover, 
specially  designed  genetic  operators  are  employed  to  make  the 
newly  generated  chromosomes  satisfy  the  constraints  automati¬ 
cally,  and  the  GA  is  accelerated  by  improving  the  initial  population 
based  on  the  optimal  control  sequence  obtained  at  the  previous 
sampling  period  and  a  local  controller.  In  addition,  the  measures  to 
achieve  offset-free  properties  are  also  discussed. 

The  paper  is  organized  as  follows:  Sections  2  and  3  introduce 
the  SOFC  system  and  some  points  that  need  to  be  considered  for 
the  design  of  a  controller.  Section  4  gives  a  brief  introduction  to  the 
SVM.  Constrained  nonlinear  predictive  controller  design  using  GA 
optimization  is  presented  in  Section  5.  The  simulations  and  discus¬ 
sions  are  presented  in  Section  6.  This  is  followed  by  the  conclusion 
in  Section  7. 

2.  SOFC  system  description 

The  SOFC  system  includes  a  fuel  processing  unit  or  reformer  and 
a  fuel  cell  stack.  Hydrogen  is  the  main  fuel  for  most  types  of  fuel 
cell.  Nevertheless,  other  fuels,  such  as  methane,  methanol,  ethanol, 
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Fig.  1.  SOFC  system  dynamic  model. 


gasoline,  and  oil  derivatives,  can  also  be  used  when  a  reformer  is 
included  in  a  fuel  cell  system  for  converting  the  fuel  into  hydrogen. 

The  basic  components  of  an  SOFC  are  the  anode,  the  cath¬ 
ode,  a  ceramic  electrolyte,  and  two  ceramic  electrodes.  In  a  fuel 
cell,  fuel  is  supplied  to  the  anode  and  air  is  supplied  to  the  cath¬ 
ode.  At  the  cathode-electrolyte  interface,  oxygen  molecules  accept 
electrons  from  the  external  circuit  to  form  oxide  ions.  The  elec¬ 
trolyte  layer  allows  only  oxide  ions  to  pass  through  it,  and  at  the 
anode-electrolyte  interface,  hydrogen  molecules  present  in  the 
fuel  react  with  oxide  ions  to  form  steam  with  the  release  of  elec¬ 
trons.  These  electrons  pass  through  the  external  circuit  and  reach 
the  cathode-electrolyte  layer,  and  thus  the  circuit  is  closed. 

In  this  paper,  a  widely  accepted  dynamic  model  of  an  SOFC  sys¬ 
tem  is  adopted  [2-7],  as  shown  in  Fig.  1,  where  Vdc  denotes  the 
stack  output  voltage  (V),  qf  the  natural  gas  flow  rate  (mol  s-1 ),  and 
/  the  external  current  load  (A);  pH2>  Po2>  and  Ph2o  denote  the  par¬ 
tial  pressures  of  hydrogen,  oxygen,  and  water  (Pa),  respectively; 
and  and  q1^  are  the  input  flow  rates  of  hydrogen  and  oxygen 
(mol  s-1 ),  respectively.  Table  1  contains  the  parameters  of  the  SOFC 
model  [2]. 

Applying  Nernst’s  equation  and  taking  into  account  ohmic,  con¬ 
centration,  and  activation  losses  (i.e.,  pohrnic,  Vconc ,  and  qact),  the 
stack  output  voltage  is  represented  as  follows  by  applying  Laplace 
transforms  [2-7]: 


V0=N0 

where 
Ph2  = 

Po2 
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1/2 
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3.  Design  considerations  for  an  SOFC  system 

The  following  aspects  should  be  considered  in  the  design  of  a 
controller  for  the  SOFC  system: 


Vdc  ~Vq  —  rjact  —  q ohmic  ~  qconc 


(1) 


Table  1 

Parameters  in  the  SOFC  system. 


Parameter  Value  Unit  Representation 


To 

1273 

K 

Fo 

96,485 

Cmol-1 

Ro 

8.314 

J  mol-1  K_1 

Eo 

1.18 

V 

N0 

384 

- 

Kr 

0.996  x  10“3 

mol  s_1  A-1 

I<H2 

8.32  x  10-6 

mol  s_1  Pa-1 

Kh20 

2.77  x  10“6 

mol  s-1  Pa- 

I<o2 

2.49  x  10“5 

mol  s-1  Pa-1 

Th2 

26.1 

s 

TH20 

78.3 

s 

X°2 

2.91 

s 

Th-o 

1.145 

- 

r 

0.126 

Q 

rf 

5 

s 

d 

0.05 

- 

P 

0.11 

- 

h 

800 

A 

Absolute  temperature 

Faraday’s  constant 

Universal  gas  constant 

Ideal  standard  potential 

Number  of  cells  in  series  in  the  stack 

Constant,  /<'r=N0/4F0 

Valve  molar  constant  for  hydrogen 

Valve  molar  constant  for  water 

Valve  molar  constant  for  oxygen 

Response  time  of  hydrogen  flow 

Response  time  of  water  flow 

Response  time  of  oxygen  flow 

Ratio  of  hydrogen  to  oxygen 

Ohmic  loss 

Time  constant  of  the  fuel  processor 
Tafel  constant 
Tafel  slope 

Limiting  current  density 


(1)  A  nonlinear  controller  is  preferred  over  a  linear  controller 
because  of  the  nonlinearity  of  the  SOFC  in  response  to  a  change 
of  operating  points.  The  stationary  voltage/current  character¬ 
istics  of  an  SOFC  system  are  depicted  in  Fig.  2,  which  shows 
that  the  SOFC  exhibits  nonlinear  behavior  over  a  wide  operat¬ 
ing  regime.  The  stack  voltage  usually  shows  significant  changes, 
especially  at  low  and  high  current  loads,  and  an  overloaded 
current  leads  to  a  rapid  deterioration  of  the  operating  stack 
voltage. 

(2)  As  a  measurable  disturbance,  the  current  load  I  should  be  uti¬ 
lized  to  construct  a  feedforward  loop  to  keep  the  output  voltage 
steady  by  speeding  up  the  initial  response  of  fuel  flow  to  drastic 
current  changes. 

(3)  Besides  ordinary  actuator  constraints  on  the  control  signal,  the 
fuel  utilization  of  the  SOFC  system  should  also  be  kept  within  a 
safe  range  for  as  long  as  possible.  Fuel  utilization  is  one  of  the 
most  important  operating  variables  that  can  affect  the  perfor¬ 
mance  of  an  SOFC.  Fuel  utilization  is  defined  according  to  Eq. 
(9): 


<-<?h2  =Qh2  _  2Krl 

<  <  < 


(9) 
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Fig.  2.  Voltage-current  characteristics  of  an  open-loop  SOFC  under  different  fuel 
flow  conditions. 


Fig.  3.  Block  diagram  of  the  proposed  nonlinear  predictive  controller. 


a  quadratic  programming  (QP)  problem  as  follows: 

m  mm 

max ^y, (at  -  ap  -  ~  ai){aI  ~  aj PK(-Xi’  xi ) 

1=1  1=1  j=  1 

m 

-e^(a++ap  (13) 

i=l 


where  q^,  q°2,  and  qj^  are  the  hydrogen  input,  output,  and 
reacted  flow  rates,  respectively.  The  desired  range  of  fuel  uti¬ 
lization  is  from  0.7  to  0.9.  This  is  in  order  to  prevent  overused 
and  underused  fuel  conditions;  an  overused-fuel  condition  can 
lead  to  permanent  damage  to  the  cells  due  to  fuel  starvation, 
while  an  underused-fuel  situation  results  in  low  efficiency  of 
the  SOFC  [3-6]. 


4.  Support  vector  machine 

This  section  briefly  reviews  the  s- SVM  algorithm  [16,17]. 

Consider  the  training  data  set  D  =  {*2,y2}™r  where  *2  is  the  ith 
input  data  in  the  input  space  Rn  and  y2  e  R  is  the  corresponding 
output  value.  The  SVM  approximates  the  relationship  between  the 
input  and  output  data  points  in  the  following  form: 


F{Xj)  =  (w,  </>(x,))  +  b  (10) 

where  w  is  a  vector  in  the  feature  space  F  c  Rn ,  0(x2 )  is  a  mapping 
from  the  input  space  to  the  feature  space  F,  b  is  the  bias  term,  and 
(., .)  denotes  the  inner  product  in  F. 

The  £-SVM  model  is  aimed  at  minimizing  the  loss  func¬ 
tion  i||w||2  +  |y2-  -  F(x2-)|g  based  on  the  following 

^-insensitive  model: 

|yi  -  F(Xi ) |g  =  max|0,  |y,-  -  F(x,- ) |  -  e]  (11) 

The  optimization  problem  can  be  formulated  as: 

m 

min  J(w,b,|+,|r)=  min  l||w||2  +  FVVf+  +  £r)  (12) 

W,b,£.  Z  in  Z - ' 

,S1  i=1 


S.t.  , 


£)(«+- ap  =  0. 


1=1 


otj~,ai  e 


0,  - 

m 


i  =  1 , . . . ,  m 


where  /C(x2,  Xj)  =  c|)r(x2)c|)(Xj)  is  a  kernel  function.  Motivated  by  Mer¬ 
cer’s  condition,  the  kernel  function  handles  the  inner  product  in  the 
feature  space  and  hence  the  explicit  form  of  0(x2  )  does  not  need  to 
be  known. 

The  solution  of  the  QP  problem,  Eq.  (13),  gives  the  optimum  val¬ 
ues  of  at  and  or.  One  then  obtains  w  =  (at  -  ar)0(x2),  where 
some  of  the  coefficients  (at  -  ap  have  non-zero  values  and  the 
corresponding  training  data  points  have  an  approximation  error 
equal  to  or  larger  than  s.  When  only  the  support  vectors  corre¬ 
sponding  to  non-zero  values  of  /32  =  (at  -  ar)  are  considered,  Eq. 
(10)  becomes: 


/ 

F(xi)  =  Xj)  +  b  (14) 

J=1 


where  l  denotes  the  number  of  support  vectors  in  the  model.  The 
bias  b  is  calculated  as: 


m 

b  =  ~l^2  (<*i  -«r)  {K(Xr,Xi)  +  K{xs,Xi))  (15) 

1=1 


where  xr  and  xs  are  support  vectors. 


5.  Constrained  model  predictive  control  with  a  terminal 
cost  based  on  genetic  optimization 


( y,-  -  (w,  0(*,O)  -  b  <  s  + 
s.t.  ^  (w, 0(Xj ))  +b-yi<s  +  $r 

U+’^° 

where  s  is  the  maximum  tolerable  error,  £+  and  are  slack 
variables,  II  •  II  is  the  Euclidean  norm,  and  C  is  a  parameter  that 
represents  a  trade-off  between  the  model  complexity  and  the  tol¬ 
erance  to  an  error  larger  than  s.  The  dual  form  of  Eq.  (12)  becomes 


The  structure  of  the  proposed  nonlinear  predictive  control  is 
given  in  Fig.  3,  where  Vr  is  the  set-point  for  the  output  voltage.  A 
nonlinear  SVM  model  is  first  identified  to  approximate  the  behavior 
of  the  SOFC  system,  and  then  a  specially  designed  GA  is  employed  to 
solve  the  resulting  nonlinear  constrained  predictive  control  prob¬ 
lem.  In  addition,  a  dynamic  constraint  unit  is  adopted  in  Fig.  3 
in  order  to  meet  the  requirements  for  fuel  utilization  and  control 
constraints. 
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5A.  Identification  of  the  SOFC  system  based  on  SVM 

To  achieve  nonlinear  predictive  control,  a  nonlinear  SVM  model 
is  used  to  approximate  the  behavior  of  the  SOFC  system. 

Considering  the  effect  of  the  dynamics  of  the  SOFC  system  up 
to  third  order,  the  output  of  the  system  at  sampling  instance  k  + 1 
can  be  reasonably  written  as  a  function  of  past  outputs,  control 
inputs,  and  measurable  disturbance  inputs  at  the  kth,  (k  -  1  )th,  and 
(/<  -  2)th  instants.  Thus: 

Vdc(k  +  1 )  =  F[ Vdc(k),  Vdc(k  -  1 ),  Vdc(k  -  2),  qf(k),  qf(k  -  1 ), 

Qf(k  -  2),  I(k),  I(k  —  1),  I(k  -  2)] 

i 

=  x(k))  +  b  (16) 

1=1 


quadratic  (LQ)  controller  as  the  local  controller  and  a  corresponding 
terminal  cost  has  been  incorporated  into  the  performance  index  of 
the  standard  GA-based  MPC. 

The  performance  index  and  the  constraints  of  the  MPC  for  the 
SOFC  system  are  defined  by: 

JV-l 

J{k)  =  ^2  [(#(k+i+l \k)-x)T Q(x{k  +  i  + 1 1 k)  -  x) 

i=0 

+  RAqj(k  +  i\k)\  +  V(x(k  +  N +  l\k)-x)  (17) 

/  x{k  +  i  +  1 1 k)  =  F'[x(k  +  i\k),  qf{k  +  i\k),  Iv{k  +  i)] 

\Vdc(k)  =  l  1  0  0  0  0 ]x(k) 


where  x{k)  =  [Vdc{k),  Vdc(k  -  1 ),  Vdc(k  -  2),  qf(k),  qf(k  -  1 ),  qf 
(/<- 2), /(/<),/(/<- 1  ),/(/<- 2)]t  is  the  input  vector  to  the  SVM 
model  at  instant  k,  Xj(i  =  l,  2,  . . .,  1)  are  support  vectors  in  the 
model,  and  ft  and  b  are  constant  coefficients  obtained  through 
training. 

The  procedure  for  SVM  identification  of  the  SOFC  system  can  be 
summarized  as  follows: 


(1)  A  set  of  measured  I/O  data  (qj.,/1,  VJC),  i  =  1, 2, . . . ,  m  is  first 
acquired  from  the  SOFC  system.  The  samples  should  cover 
the  whole  working  range.To  prevent  some  elements  that  have 
larger  original  absolute  values  from  dominating  the  final  kernel 
value,  it  is  necessary  to  carry  out  some  preprocessing  of  the  raw 
data  before  feeding  them  into  the  SVM  model.  In  this  project,  all 
of  the  feature  elements  and  the  target  values  have  been  scaled 
so  that  they  fall  into  the  range  of  [-1,1  ].  When  using  this  SVM 
model,  the  computed  target  value  should  be  converted  back  to 
the  same  scales  that  were  used  for  the  original  target  values. 

(2)  Next,  the  training  samples  of  the  SVM  are  organized  using  xt  = 

r  n  T 


ji+ 2  ji+ 1  ji 


as  the  input  vari¬ 


able  andyz  =  Vj+3  as  the  output  variable,  i  =  1,  2, . . .,  (m  -3). 

(3)  The  parameters  of  the  SVM,  including  the  weighting  coeffi¬ 
cient  C,  the  width  parameter  or  (a  Gaussian  kernel  function  is 
adopted),  and  s,  are  set.  a  and  C  affect  the  generalization  ability 
of  the  SVM  in  opposite  directions.  Setting  C  too  low  will  result 
in  insufficient  learning  from  the  training  data,  while  setting  C 
too  high  will  lead  to  overfitting. 

(4)  Lastly,  the  optimization  problem,  Eq.  (12),  is  constructed  and 
solved. 


Flerein,  the  optimization  problem,  Eq.  (12),  is  practically  solved 
by  its  dual  problem,  Eq.  (13),  using  the  sequential  minimal  opti¬ 
mization  (SMO)  algorithm  [18,19].  The  values  of  at  and  af  are 
then  determined,  and  the  nonlinear  SVM  model  of  the  SOFC  can  be 
acquired  according  to  Eq.  (16). 

5.2.  Formulation  of  the  predictive  control  problem  for  an  SOFC 

Because  of  the  use  of  a  nonlinear  SVM  predictive  model,  an  NP- 
hard  nonlinear  optimization  problem  needs  to  be  solved  to  achieve 
the  MPC.  A  GA  has  been  shown  to  give  better  performance  than 
a  quasi-Newtonian  optimization  technique  in  solving  this  kind  of 
optimization  problem  [12].  In  this  study,  the  standard  GA-based 
model  predictive  control  has  been  modified  by  taking  advantage 
of  the  well-developed  theory  on  the  stability  of  model  predictive 
control  [20-22]. 

The  three  key  ingredients  of  a  stabilizing  predictive  control  are 
summarized  in  Ref.  [20],  which  comprise  a  terminal  set,  a  terminal 
cost,  and  a  local  controller.  In  this  study,  we  have  employed  a  linear 


Qf  min  —  Qf(k  +  i\k)  <  C[f  max?  i  —  0,  .  .  .  ,  N  1  (19) 

Acfy-  min  —  Aqftk  +  i\k)  <  Aqjmax,  i  =  0, . . . ,  N  —  1  (20) 

where  N  is  the  prediction  horizon;  F( . )  is  the  state-space  represen¬ 
tation  of  the  SVM  model  F(.)  in  Eq.  (16)  by  setting  the  state  vector 
*(/<)  =  [Vdc(k\  Vdc(k  -  1 ),  Vdc(k  -  2),  qfik  -  2),  qfik  -  1  )]T ;  x(k  +  i  +  l|lc) 
is  the  predicted  state  at  instant  k  +  i  + 1  based  on  the  current  state 
x{k)  and  control  sequence;  Aqfik  +  i\k)  =  qj(k  +  i\k)  -qfik  +  i-  1 1/<)  is 
the  predicted  change  in  the  control  input  signal;  x  is  the  equi¬ 
librium  value  of  the  state  vector  corresponding  to  the  current 
set-point  (see  (27));  Q=Qr >0,  R> 0;  [q/min,  q/max]  are  dynamic 
constraints  to  meet  the  requirements  for  fuel  utilization  and  con¬ 
trol  constraints,  which  will  be  described  in  Section  5.3;  and  Iv(k)  = 
[/(k),/(k-l),/(/c-2)]r. 

The  last  term  in  Eq.  (17)  represents  the  terminal  cost.  This  ter¬ 
minal  cost  represents  the  stabilizing  cost  that  would  be  required 
when  the  system  is  to  be  controlled  beyond  the  finite-time  horizon 
N  toward  the  infinite-time  horizon.  We  assume  that  this  stabilizing 
controller  can  be  designed  by  a  fictitious  local  LQcontroller  near  the 
equilibrium  point. 

The  optimal  cost  in  driving  the  system  to  equilibrium  from  the 
time  instant  k  +  N+ 1  to  the  infinite  time  is  defined  by: 

^(x(k  +  N  +  1  \k)  -  x)  =  (x(k  +  N  +  1  \k)  -  x)TP(x(k  +  N  +  1  |k)  -  x) 

(21) 

where  P  is  the  symmetric  positive  semi-definite  solution  of  the 
algebraic  Riccati  equation: 

P  =  AtPA-AtPB(BtPB  +  rCbtPA  +  Q  (22) 

where  A  and  B  are  constant  matrices  of  appropriate  dimensions 
obtained  by  linearizing  or  identifying  the  SOFC  system  Eq.  (1) 
around  the  nominal  operation  point;  and  Q.  =  Q.  >0  and  R  >  0 
are  the  weighting  matrices  for  the  local  LQcontroller. 

The  local  controller  JCL q  corresponding  to  the  terminal  cost  is 
designed  as: 

tfLQ  =  -( BtPB  +  R)~'btPA  (23) 

In  formulating  the  optimization  problem,  (17)-(20),  we  assume 
that  there  is  a  finite  horizon  length  N,  such  that  the  prediction  of  the 
state  vector,  x(k  +  N+ 1 1 k)  e  £2,  where  £2  is  the  terminal  set,  and  the 
constraints  can  be  assumed  to  be  inactive  for  i  >  N,  because  it  is  very 
difficult  to  calculate  the  terminal  set  for  a  given  nonlinear  system, 
and  also  to  reduce  the  computational  load.  As  a  result,  the  terminal 
set  constraint  is  omitted  in  (17)-(20).  Note  that  the  assumption 
concerning  the  terminal  set  may  be  met  if  the  prediction  horizon 
is  long  enough  or,  in  the  case  of  short  horizon,  a  large  weighting 
matrix  R  is  used  in  designing  the  local  controller. 

The  local  LQ  controller  is  only  used  to  determine  a  terminal 
penalty  matrix  P  off-line  and  to  help  the  GA  improve  the  initial 
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population,  which  will  be  illustrated  in  Section  5.5.  This  method 
does  not  require  the  globally  optimal  input  profile  to  be  found 
numerically  at  every  step.  Stability  only  requires  feasible  solutions 
to  the  optimization  problem.  The  computational  (and  performance) 
advantage  of  this  scheme  lies  in  the  fact  that  shorter  horizons  can  be 
used,  without  jeopardizing  performance  and  stability.  This  is  espe¬ 
cially  beneficial  for  on-line  application  of  the  GA-based  predictive 
control  approach. 

5.3.  Calculation  of  the  dynamic  constraints 

We  assume  that  the  safe  range  of  fuel  utilization  is  within  [Ufm[n, 
u/maxL  and  that  the  actuator  constraints  on  the  control  signal  are 
[9/min*  q/max ]•  The  dynamic  constraint  unit  in  Fig.  3  is  designed  as 
follows. 

In  terms  of  Eq.  (9),  to  guarantee  a  safe  utilization,  the  hydro- 
gen  flow  should  be  within  [qH2min,  <?H2max].  where  qH2min(k)  = 
2Kr/(k)/uf 

max  and  (}h2  max(k)  =  2 KrI(k)/uf  min- 

Because  it  is  appropriate  to  approximate  the  fuel¬ 
processing  unit  using  a  first-order  transfer  function,  we 
can  calculate  the  constraints  [q'f  .  ,q'  1  on  the  fuel 

flow  (the  natural  gas  flow  rate)  from  [gH2min>  4H2max]-  For 
example,  in  the  case  of  sample  time  Ts  =  ls,  QfmaxU<)  = 
(<7H2max(k)  -  0.8187(jH2(k))/0.1813,  q’  (k)  =  (<Jh2  min(k)  - 
0.8187<jH2(k))/0.1813.  ' 

The  final  dynamic  constraints  [q/min,  Qf  maxi on  the  fuel  flow  can 
thus  be  calculated  using: 

qf  min  =  max(q/min,  q’fmin)  (24) 

Qf  max  —  min(qfmax>  4/ max)  (25) 

5.4.  Offset- free  output  tracking 

To  achieve  offset-free  output  tracking,  a  simple  disturbance  esti¬ 
mator  is  run  at  each  step,  assuming  that  the  state  disturbance  is 
given  by  the  difference  between  the  latest  measured  state  and  the 
previously  expected  state: 

d(/<)  =  x{k)  -  x(k\k  -  1)  (26) 

This  is  then  used  to  estimate  the  equilibrium  values  of  the  state 
vector  x  and  the  fuel  flow  qf,  assuming  that  the  disturbance  will 
remain  constant  at  this  estimated  value:  d  =  d{k).  Specifically,  x 
and  qf  are  found  as  the  solutions  to  the  simultaneous  equations: 


5.5.  Genetic  optimization  of  control  inputs 

Due  to  the  use  of  an  SVM  model,  the  proposed  method  for¬ 
mulates  a  dynamic  nonlinear  optimization  problem,  to  which  the 
conventional  optimization  techniques  cannot  be  easily  applied. 
Therefore,  in  this  work,  the  on-line  optimization  problem  is  solved 
using  a  GA. 

To  deal  with  constraints  in  this  optimization  problem,  the 
penalty  function  method  is  commonly  used.  However,  this 
approach  lowers  the  efficiency  of  a  GA,  because  of  the  waste  of 
genetic  material  due  to  unfeasible  solutions  stemming  from  stan¬ 
dard  genetic  operators  in  the  population.  In  this  study,  specially 
designed  genetic  operators  have  been  employed  to  make  the  newly 
generated  chromosomes  satisfy  the  constraints  automatically. 

The  algorithm  starts  with  an  initial  population  of  chromosomes, 
which  represent  possible  solutions  of  the  optimization  problem. 
The  objective  function  is  computed  for  each  chromosome.  New 
generations  are  produced  by  the  genetic  operators,  which  are  des¬ 
ignated  as  selection,  crossover,  and  mutation.  The  algorithm  stops 
after  the  maximum  permissible  time  has  elapsed. 

A  chromosome  that  is  a  candidate  solution  of  the  optimiza¬ 
tion  problem  is  represented  by  s*,  the  elements  of  which  comprise 
present  and  future  control  inputs,  and  has  the  following  structure 
[13]: 

Si  =  lqfi(k)  qfi(k  + 1)  •••  qfi(k  +  N- 1)],  i  =  (29) 

where  k  indicates  the  current  time,  and  L  is  the  number  of  chromo¬ 
somes.  The  algorithm  can  be  described  as  follows: 

Step  1  (initial  population  generation):  the  number  of  iterations 
iter=\  is  set.  Predictive  control  uses  the  receding  horizon  princi¬ 
ple,  which  implies  that  an  evolution  has  to  be  calculated  at  each 
time  step.  Hence,  the  past  evolutions  give  important  information 
that  can  be  used  to  improve  the  initial  population  of  the  current 
evolution 

The  optimal  input  sequence  obtained  at  k  - 1  is  assumed  to 
be  U*(k  -  1)  =  [q*(k  -  1),  q*(k),  ...,  q*(k  +  N  -  2)}.  At  the  cur¬ 

rent  time  k,  we  consider  a  “shifted”  input  sequence  £/(/<)  as  shown 
below,  where  the  last  gene  takes  qf(k  +  N  -  1),  obtained  using  the 
designed  local  controller,  to  be  one  of  the  initial  chromosomes.  This 
chromosome  might  be  a  very  good  guess  for  the  solution  of  the  next 
optimization  problem. 

=  qf(k ),  ...,  qf(k  +  N-3),  qf(k  +  N-2)\ 

l  applied  J 


( x  =  F'[x,qf,Iv(k)]  +  d  r27s 

|Vr  =  [l  0  0  0  0]x  1  ; 

Note  that  this  results  in  offset-free  control  in  the  presence  of  an 
unknown  but  constant  disturbance,  even  if  the  steady-state  gains 
in  the  model  are  not  accurate. 

Another  measure  to  achieve  offset-free  properties  is  to  intro¬ 
duce  an  integration  action.  We  assume  that  x'{k)  =  [AVdc(k), 
A Vdc{k-\),  AV{k-2),  Aqj(k-2)t  Aqf(k- 1),  e{k)]T,  where 
e(/<)  =  Vr(k)  -  Vdc(k).  The  performance  index  is  accordingly  modified 
as: 

N-i 

m  =  J2lx'T<~k  + '  +  1  lk)Q*'(k  +  i  +  1  Ik)  +  RAqj(k  +  i\k)] 

1=0 

+  vI'(x/(k  +  N  +  l|/<))  (28) 

Because  this  method  may  result  in  integrator  wind-up,  this  part 
of  the  contents  is  omitted  herein. 


*,=(?(*)  =  Ly(*),  qfk  + 1),  ....  f(k  +  N-2),  qf(k  +  N-\) 

L  new  tail 

where  the  newly  added  tail  is  defined  by 
(  low  if  q'j  <  low 


qf(k  +  N-  1)=<  Q}  if  low  <  q'f  <  upper 

(30) 

^  upper  if  q'j:  >  upper 

q'f  =  Kiofx{k  +  N  -  1 1/<)  -  x)  +  qf 

(31) 

low  =  max{0,  qj(1<  +  N  -  2)  +  Aq/min} 

(32) 

upper  =  min{l ,  q*[k  +  N  -  2)  +  Aq/max} 

(33) 

To  satisfy  both  the  control  and  control  move  constraints,  we  use 
a  simple  procedure  to  generate  the  remaining  L  -  1  chromosomes 
s2  -  sL  of  the  initial  population. 
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(1 )  For  the  current  input  value: 

Qfi(k)  =  min(q/max,  max(q/min,  rcmd[qf(k  -  1) 

Aq^  min’  Qf(k  —  1 )  +  Aq/max]))’  2  <  i  <  L 


Fors"™,  it  is  supposed  that  d  =  qfl(k+z)-q^,+i)(k+z- 1),  and 
then: 


qfl(.k+Z+j)= 


qfi(k+z+j)-(d-Aqf 

max  ),  if  d  >  A qf 

max 

<lfl(.k+z+j)-{d- Aqf 

min  ),  if  d  <  Aqfmin 


(2)  For  the  rest  of  the  input  values: 


(0<j<N  —  z— 1) 


(36) 


QfiU<  +J)  =  min(g/max,  max(q/min,  rand[qfi{k  +j  -  1) 

+  min?  + J  —  1)  +  max  D), 

2  <i  <L,  1  <  j  <  N  —  1 


In  the  above  equations,  rand  [.]  is  a  random  number  within  [.], 
and  a  new  random  number  is  generated  each  time. 

Step  2  (fitness  function  evaluation):  the  objective  function  of  Eq. 
(17)  is  evaluated  for  all  of  the  chosen  chromosomes.  Their  fitness 
value  is  then  calculated  according  to: 

fitness(i )  =  >  i  =  1 , 2, . . . ,  L  (34) 

where  fi  is  the  value  of  the  objective  function  for  the  ith  chromo¬ 
some. 

The  normalized  fitness  value  of  each  chromosome,  that  is,  the 
selection  probability,  is  then  calculated  by: 

p,=  fneSS(0  (35) 

Step  3  (selection  operation):  the  m2  best  individuals  are 
retained  and  reintroduced  into  the  population  of  the  next  gen¬ 
eration.  Therefore,  the  partly  optimized  chromosomes  will  not 
get  lost  in  spite  of  the  disruption  of  building  blocks  during 
crossover. 

The  rest  of  the  L-m2  chromosomes  are  generated  according  to 
their  selection  probabilities.  The  chromosomes  with  high  fitness 
value  have  a  greater  chance  of  being  selected. 

Step  4  (crossover  operation):  for  each  chromosome,  a  random 
number  r\  between  0  and  1  is  generated.  If  r\  is  lower  than  the  prob¬ 
ability  of  crossover  pc,  this  particular  chromosome  will  undergo 
the  process  of  crossover,  otherwise  it  will  remain  unchanged.  The 
selected  chromosomes  are  paired  and  for  each  selected  pair  one  of 
the  following  two  crossover  operations  is  implemented  with  equal 
probability 

(1 )  The  one-point  crossover  operation 

A  random  integer  z  between  1  and  N-  1  is  generated,  which 
indicates  the  position  of  the  crossing  point.  Two  new  chromo¬ 
somes  are  produced  by  interchanging  all  of  the  members  of  the 
parents  following  the  crossing  point,  which  can  be  expressed 
graphically  as  follows: 


Similar  equations  can  be  obtained  for  the  chromosome  s 7ew. 
(2)  The  uniform  crossover  operation 

For  the  uniform  crossover  operation,  two  new  chromosomes 
based  on  s,  and  sI+1  are  produced  by: 


/  s"ew  =  r2  •  st  +  (1  -  r2)  •  sI+1 
\s['+7  =  (l  -  r2)  ■  Sj  +  r2  ■  s;+1 

where  r2  is  a  random  number  between  0  and  1 . 


(37) 


Step  5  (mutation  operation):  for  every  member  of  each  chromo¬ 
some,  a  random  number  r3  between  0  and  1  is  generated.  If  r3  is 
lower  than  the  probability  of  mutation  pm,  this  particular  member 
of  the  chromosome  will  undergo  the  process  of  mutation;  other¬ 
wise,  it  will  remain  unchanged.  For  the  selected  members,  lower 
and  upper  bounds  [bi(j),bu(j)]  are  defined  as: 

r  min(A qfmax+qf(k  -  1),  Aqfmax+qfi(k+l),  q/max),  j  =  0 
bu(j)=<  mm(Aqfmax+qfi(k+j-\),  Aqfmax+qfi(k+j+\),qfmax),  0<j<N-l  (38) 
l  min(Aq/max+^(k+j-l),q/max),  j=N- 1 

f  max(Aqfmin+qf(k~L),  Aqfmin+qfi(M),qfmin),  j  =  0 
bi(j)=  <  max( Aqfmin+qfi(k+j-l),  Aqfmin+qfl(k+j+l),  qfmi n),  0  <  j  <  N  -  1  (39) 

^  max(A qfmin+qfi(k  +j  -  1),  q/min),  j  =  N  -  1 

The  above  bounds  define  the  region  of  values  that  will  produce 
a  feasible  solution.  The  mutation  operation  is  then  achieved  by  the 
generation  of  a  random  number  within  [b/(j),bu(j)]. 

q™w(k+j)  =  rm(j)  if  r3<pm  (40) 

where  rm(j)  is  a  random  number  within  [bi(j),bu(j)]- 

Step  6  (repeat  or  stop):  if  the  maximum  allowed  time  has  not 
elapsed,  the  algorithm  is  set  and  returned  to  Step  2.  Otherwise,  the 
algorithm  is  stopped  and  the  chromosome  that  produced  the  high¬ 
est  value  of  the  fitness  function  throughout  the  entire  procedure  is 
selected. 

The  previous  modified  GA  makes  it  possible  to  calculate  sub- 
optimal  control  in  real  time. 


6.  Simulation  results 


In  this  section,  we  describe  application  of  the  proposed  nonlin¬ 
ear  predictive  controller  to  the  SOFC  problem  to  achieve  safe  fuel 
utilization  and  maintain  operational  constraints  when  only  the 
voltage  output  is  measurable  on-line.  The  sampling  rate  of  the  SOFC 
and  the  MPC  was  chosen  as  Ts  =  1  s.  For  all  simulations,  the  following 
parameters  were  set:  [q/min,  q/max]  =  [0, 1.2]  mol  s~\  A q/max  = 
-Aq/min  =  0.7  mol  s-1,  [u/min,  u/max]  =  [0.7,  0.9],  pc  = 

0.8,  Pm  —  0.1,  L  =  20,  m2  =2,  ft  =  10,  Q  =  Q  = 
diag{  0.1  0  0  0  0). 


Si  =  {qfi(k), . . . ,  qfi(k  +  z  -  1 ),  \qfi(k  +  z), . . . ,  qfi(k  +  N-  1 )} 
s,-+1  =  {(Z/(l+i )(/<),  •  •  • ,  <Z/(i+i)(fc  +  z-  1 ),  |<2/(i+i)(fc  +  z), . . . ,  q/(l+i)(/<  +  N-  1 )} 
the  one  point  crossover 

s- ew  =  {qfi(k), . . . ,  qfi(k  +  z  -  1 ),  \qm^(k  +  z), . . . ,  qf{i+^(k  +  N-  1 )} 

=  (<7/(i+ nOO,  •  •  • ,  <?/(i+i)(k  +  z  -  1 ),  | qfi(k  +  z), . . . ,  qfi(k  +  N-  1 )} 

The  previous  operation  might  produce  infeasible  offspring  if 
the  input  values  at  the  crossing  point  do  not  satisfy  the  control 
move  constraints.  This  situation  is  avoided  by  the  following  cor¬ 
rection  mechanism  for  each  of  the  new  chromosomes  s"ew  and 
snew,  which  modifies  the  values  of  the  input  parameters  after 
the  crossing  position  so  that  the  control  move  constraints  are 
satisfied. 


6.1.  Identification  of  the  SOFC  system 

Open-loop  input-output  data  are  required  to  train  the 
SVM  model  of  an  SOFC  system.  Open-loop  input-output  data 
samples  were  obtained  by  exciting  the  open-loop  SOFC  sys¬ 
tem  using  designed  sine  signals  0.7823 +  0.3  sin(0.5t)  sin  t  and 
300  +  50sin(0.03t)sin(0.04t)  for  the  fuel  and  the  current  demand 
(i.e.,  qf  and  /),  respectively,  which  were  collected  over  1 000  s  and  are 
plotted  in  Fig.  4.  All  1000  of  the  sampled  data  points  were  divided 
into  two  sets,  the  training  set  and  the  testing  set,  where  the  train¬ 
ing  set  contained  900  data  points  and  the  testing  set  contained  the 
remaining  100  data  points. 
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Fig.  4.  Input  excitation  and  output  response  signals  of  the  SOFC. 
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Fig.  5.  Training  and  test  results  of  the  SVM  model. 


We  first  identified  the  SVM  model  of  the  SOFC  system  with  dif¬ 
ferent  parameter  settings,  giving  the  results  shown  in  Table  2.  The 
quality  of  the  approximation  was  assessed  using  the  root-mean- 
square  error  (RMSE)  between  the  samples  and  the  SVM  model. 

In  order  to  make  a  trade-off  between  the  accuracy  and  the  com¬ 
plexity  of  the  SVM,  we  chose  C  =  1 0,000,  s  =  0.01 ,  and  a  =  40  in  terms 
of  the  test  error  and  the  number  of  support  vectors  in  the  following 
simulations.  The  corresponding  training  and  test  results  are  plotted 
in  Fig.  5,  which  shows  that  the  SVM  can  approximate  the  behavior 
of  the  SOFC  system  with  good  accuracy. 


6.2.  Predictive  control  of  the  SOFC  system 


By  linearizing  Eq.  (1 )  around  its  equilibrium  point,  we  obtained 
the  transfer  function  of  the  SOFC  system: 

G(s)=%(£)=k  - — — — — TT7 - ^  (41) 

q/(s)  g  (TfS  +  l)(tH2S  +  1)(To2S  +  1) 

where  kg  =  Y1+Y2,  r  =  (Y,t02  +  Y2Th2'I/(Y^  +  Y2),  Y,  = 

NoRoW(2FoJ4i2PH2,o),  ^2  =  NoRo^o/(4ForH-o^o2Po2,o)  • 

Ph2,o  and  Po2,o  denote  the  steady-state  partial  pressures  of 
hydrogen  and  oxygen,  respectively. 

Specifically,  around  the  nominal  operation  point  with 
qf=  0.7023  mol  s-1 ,  /  =  300  A,  Vdc  =  333  V: 


G(s)  =  230.37  • 


5.85s +  1 

(5s  +  1)(26.1s  +  1)(2.91s  +  1) 


(42) 


Table  2 

Training  and  test  results  with  different  parameter  settings. 


o 

e 

C 

Training  error 
(RMSE) 

Test  error 
(RMSE) 

Number  of 
support  vectors 

20 

0.01 

10,000 

0.0462 

0.0454 

24 

40 

0.01 

10,000 

0.0435 

0.0391 

16 

60 

0.01 

10,000 

0.0416 

0.0411 

14 

40 

0.02 

10,000 

0.0985 

0.0888 

14 

40 

0.01 

1000 

0.0458 

0.0405 

19 

Converting  Eq.  (42)  to  its  equivalent 

discrete  state-space  form, 

we  obtain: 

J  x(k  +  1 )  =  Ax{k)  +  BAqf(k) 

(43) 

\  Vdc(k)  =  Cx(k) 

"2.49 

-2.051 

0.5588 

-1.141 

1.5981" 

1 

0 

0 

0 

0 

where  A  - 

0 

1 

0 

0 

0 

, 

0 

0 

0 

0 

1 

_  0 

0 

0 

0 

1 

B  =  [1.553  0 

0  0 

HT, 

C  =  {  1 

0  0  0 

0] 

The  following  positive  semi-definite  matrix  and  local  controller 
are  then  calculated  by  Eqs.  (22)  and  (23)  with  R  =  10,  000. 


~  88.1 

-125.9 

45.6 

-93.1 

237.1  - 

-125.9 

181 

-65.7 

134.2 

-348.6 

P  = 

45.6 

-65.7 

23.9 

-48.8 

128.6 

-93.1 

134.2 

-48.8 

99.7 

-262.6 

_  237.1 

-348.6 

128.6 

-262.6 

1219.3  _ 

I<lq  =  [  0.0318  -0.0466  0.0172  -0.0351  0.1461] 

To  illustrate  the  effectiveness  of  the  proposed  nonlinear  pre¬ 
dictive  controller,  we  assume  that  a  load  disturbance  causes  step 
changes  of  the  current  at  t  =  1 0  s  and  t = 60s,  respectively.  Let  N = 2. 
Fig.  6  shows  a  comparison  of  the  closed-loop  response  of  the  SOFC 
system  with  terminal  cost  (solid  line)  and  without  terminal  cost 
(dashed  line).  These  results  show  that  better  control  performance 
can  be  achieved  by  considering  terminal  cost,  even  though  a  short 
prediction  horizon  is  adopted.  The  main  advantage  of  adopting  a 
short  horizon  is  that  the  on-line  computational  burden  is  effec¬ 
tively  decreased  since  both  the  number  of  decision  variables  and 
the  number  of  prediction  steps  are  reduced. 

6.3.  Effect  of  the  dynamic  constraints 

To  show  the  effect  of  the  dynamic  constraints,  we  replaced 
the  dynamic  constraints  on  the  control  input  signal  in  Fig.  6  with 
the  usual  maximum  and  minimum  constraints,  which  yielded  the 
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the  dynamic  constraints.  Note  that  because  of  the  slow  response  of 
the  fuel  flow  to  the  output  voltage  and  constraints  on  the  fuel  flow, 
the  fuel  utilization  could  not  be  definitely  restricted  to  the  desired 
range  for  large  and  sudden  current  load  changes. 

7.  Conclusion 

We  have  proposed  a  nonlinear  predictive  control  strategy  to 
solve  the  SOFC  control  problem,  whereby  a  nonlinear  SVM  model 
is  first  identified  to  approximate  the  behavior  of  the  SOFC  sys¬ 
tem  using  an  SMO  algorithm,  and  then  a  specially  designed  GA 
is  employed  to  solve  the  resulting  nonlinear  constraint  predictive 
control  problem.  Meanwhile,  the  standard  performance  index  has 
been  modified  by  incorporating  a  terminal  cost,  which  makes  it  pos¬ 
sible  to  adopt  short  horizons  while  maintaining  a  satisfactory  level 
of  performance.  This  is  especially  beneficial  for  on-line  application 
of  GA-based  predictive  control.  Moreover,  the  GA  was  accelerated 
by  improving  the  initial  population  based  on  the  optimal  control 
sequence  obtained  at  the  previous  sampling  period  and  a  local 
controller.  Simulation  results  on  the  SOFC  system  have  illustrated 
that  the  proposed  method  can  successfully  deal  with  the  control 
and  control  move  constraints,  and  that  a  satisfactory  closed-loop 
performance  can  be  achieved. 


Fig.  6.  Closed-loop  response  of  the  SOFC  with  (solid  line)  and  without  terminal  cost  Acknowledgment 

(dashed  line). 
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Fig.  7.  Simulation  results  without  dynamic  constraints. 


results  shown  in  Fig.  7.  Comparing  Fig.  7  with  Fig.  6,  it  is  clear  that 
the  fuel  utilization  could  no  longer  be  kept  within  the  desired  safe 
range  for  the  current  disturbance  from  t= 62  s  to  t  =  67  s  without 


This  work  is  supported  by  the  National  Natural  Science  Founda¬ 
tion  of  China  (nos.  51076027  and  51036002). 

References 

[1]  Y.H.  Li,  S.S.  Choi,  S.  Rajakaruan,  IEEE  Trans.  Energy  Convers.  20  (2)  (2005) 
381-387. 

[2]  V.  Knyazkin,  L.  Soder,  C.  Canizares,  IEEE  PowerTech  Conference,  Bologna,  Italy, 
June  23-26,  2003,  pp.  328-333. 

[3]  X.  Wang,  B.  Huang,  T.  Chen,  J.  Process  Control  17  (2)  (2007)  103-114. 

[4]  H.B.  Huo,  X.J.  Zhu,  W.Q.  Hu,  H.Y.  Tu,  J.  Li,  J.  Yang,  J.  Power  Sources  185  (1 )  (2008) 
338-344. 

[5]  X.J.  Wu,  X.J.  Zhu,  G.Y.  Cao,  H.Y.  Tu,  J.  Power  Sources  179  (1)  (2008)  232-239. 

[6]  F.  Jurado,  J.  Power  Sources  158  (1)  (2006)  245-253. 

[7]  X.W.  Zhang,  S.H.  Chan,  H.K.  Ho,  J.  Li,  G.  Li,  Z.  Feng,  Int.  J.  Hydrogen  Energy  33 
(2008) 2355-2366. 

[8]  J.  Yang,  X.  Li,  H.G.  Mou,  L.  Jian,  J.  Power  Sources  193  (2)  (2009)  699-705. 

[9]  X.C.  Xi,  A.N.  Poo,  S.K.  Chou,  Control  Eng.  Pract.  15  (8)  (2007)  897-908. 

[10]  W.  Naeem,  R.  Sutton,  J.  Chudley,  F.R.  Dalgleish,  S.  Tetlow,  Int.  J.  Control  78  (14) 
(2005) 1076-1090. 

[11]  U.  Yuzgec,  Y.  Becerikli,  M.  Turker,  ISA  Trans.  45  (4)  (2006)  589-602. 

[12]  S.C.  Shin,  S.B.  Park,  Electron.  Lett.  34  (20)  (1998)  1980-1981. 

[13]  H.  Sarimveis,  G.  Bafas,  Fuzzy  Sets  Syst.  139  (1)  (2003)  59-80. 

[14]  C.  Onnen,  R.  Babuska,  U.  Kaymak,  J.M.  Sousa,  H.B.  Verbruggen,  R.  Isermann, 
Control  Eng.  Pract.  5  (10)  (1997)  1363-1372. 

[15]  M.G.  Na,  B.R.  Upadhyaya,  IEEE  Trans.  Nucl.  Sci.  53  (4)  (2006)  2318-2327. 

[16]  B.  Scholkopf,  A.J.  Smola,  Learning  with  Kernels:  Support  Vector  Machines,  Reg¬ 
ularization,  Optimization  and  Beyond,  MIT  Press,  Cambridge,  MA,  2002. 

[17]  A.J  Smola,  B.  Scholkopf,  A  tutorial  on  Support  Vector  Regression,  NeuroCOLT 
Technical  Report  no.  NC-TR-98-030,  Royal  Holloway  College,  University  of  Lon¬ 
don,  1998. 

[18]  G.W.  Flake,  S.  Lawrence,  Mach.  Learn.  46  (2002)  271-290. 

[19]  S.K.  Shevade,  S.S.  Keerthi,  C.  Bhattacharyya,  K.R.K.  Murthy,  IEEE  Trans.  Neural 
Networks  11  (5)  (2000)  1188-1193. 

[20]  D.Q,  Mayne,  J.B.  Rawlings,  C.V.  Rao,  P.O.M.  Scokaert,  Automatica  36  (6)  (2000) 
789-814. 

[21  ]  J.M.  Maciejowski,  Predictive  Control  with  Constraints,  Prentice-Hall,  New  York, 
2002. 

[22]  H.  Chen,  F.  Allgower,  Automatica  34  (10)  (1998)  1205-1217. 


